#!/usr/bin/perl
my $genome=$ARGV[0];
my $sample_name=$ARGV[1];  ## provide results' bam's prefix
my $file_1=$ARGV[2];
my $file_2=$ARGV[3];

my $sample = $sample_name;
my $BisSNPVersion="BisSNP-0.66.jar";
my $cpu_cores=12;
my $mem="16g";
my $genome_version=$genome;
my $readGroup_name = $sample_name;
my $dbsnp;
my $indel_1;
my $indel_2;
my $ref;
my $chr="";

if($genome_version eq "hg18"){
	$dbsnp="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/dbsnp_135.hg18.sort.vcf";
	$indel_1 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/1000G_phase1.indels.hg18.sort.vcf";
	$indel_2 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/Mills_and_1000G_gold_standard.indels.hg18.sites.sort.vcf";
	$ref="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/hg18/hg18_unmasked.plusContam.fa";
}elsif($genome_version eq "hg19"){
	$dbsnp="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/dbsnp_135.hg19.sort.vcf";
	$indel_1 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/1000G_phase1.indels.hg19.sort.vcf";
	$indel_2 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/Mills_and_1000G_gold_standard.indels.hg19.sites.sort.vcf";
	$ref="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/hg19/hg19_rCRSchrm.fa";
	
}elsif($genome_version eq "GRCh37"){
	$dbsnp="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/dbsnp_135.b37.vcf";
	$indel_1 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/1000G_phase1.indels.b37.sort.vcf";
	$indel_2 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/Mills_and_1000G_gold_standard.indels.b37.sites.sort.vcf";
	$ref="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/GRCh37/GRCh37-lite.fa";

}elsif($genome_version eq "mm9"){
	#$dbsnp="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/mouse-20111102-snps-all.annotated.mm9.vcf";
	#$indel_1 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/mouse-20110602-callable-dinox-indels.annot.mm9.vcf";
	$dbsnp="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/mouse-20111102-snps-129s1-cast.annotated.mm9.vcf";
	$indel_1 = "/export/uec-gs1/laird/shared/research/yaping/data/genome_data/dbSNP/mouse-20110602-indels-129s1-cast.annotated.mm9.vcf";

	$ref="/export/uec-gs1/laird/shared/research/yaping/data/genome_data/genome/mouse/mm9.fa";
}

my $input_bam_raw_withoutheader_super=$sample_name.".bam";

my $input_bam_raw_withoutheader=$input_bam_raw_withoutheader_super;
#if($chr ne ""){
#	$input_bam_raw_withoutheader =~ s/bam//;
#	$input_bam_raw_withoutheader .= "$chr.bam";
#}
my $input_bam_raw = $input_bam_raw_withoutheader;
$input_bam_raw =~ s/.bam//;
$input_bam_raw .= ".header.bam";
my $input_bam_realign = $input_bam_raw;
$input_bam_realign =~ s/.bam//;
$input_bam_realign .= ".realign.bam";
my $input_bam_realign_mdups = $input_bam_realign;
$input_bam_realign_mdups =~ s/.bam//;
$input_bam_realign_mdups .= ".mdups.bam";
my $input_bam_realign_mdups_recal = $input_bam_realign_mdups;
$input_bam_realign_mdups_recal =~ s/.bam//;
$input_bam_realign_mdups_recal .= ".recal.bam";

my $indel_target_interval = $input_bam_raw;
$indel_target_interval =~ s/.bam//;
$indel_target_interval .= ".indels.intervals";
my $recalFile_before=$input_bam_realign_mdups_recal;
$recalFile_before =~ s/.bam//;
$recalFile_before .= ".beforeRecal.txt";
my $recalFile_after=$input_bam_realign_mdups_recal;
$recalFile_after =~ s/.bam//;
$recalFile_after .= ".afterRecal.txt";

#my $bedfile="/home/uec-00/yapingli/software/gatk_read/Sting/resources/Normal_SNP.filterCNV.txt.minReadCov_1.snpcaller.bed";
my $tmp_dir="/export/uec-gs1/laird/shared/research/yaping/data/tmp/";
my $metrics_file=$input_bam_realign_mdups.".dupsMetrics.txt";


my $header = "#PBS -q lairdprio\n";
$header .= "#PBS -l walltime=300:00:00,mem=$mem,nodes=1:ppn=$cpu_cores\n";
$header .= "#PBS -N BisSNP-preprocess-bam-$sample-$chr\n";
#$header .= "#PBS -W depend=afterany:1550139.hpc-pbs\n";
$header .= "#PBS -j oe\n";
$header .= "cd \$PBS_O_WORKDIR\n";
#merge bam
#print "$ARGV[0]\t$ARGV[1]\t$ARGV[2]\t$ARGV[3]\n";
if($#ARGV >= 3){
	$header .= "samtools merge -f $input_bam_raw_withoutheader_super $file_1 $file_2\n";
	$header .= "samtools index $input_bam_raw_withoutheader_super\n";
}


#reorder bam
#$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/picard/ReorderSam.jar ";
#$header .= "I=$input_bam_raw_withoutheader_super_no_order O=$input_bam_raw_withoutheader_super R=$ref S=true CREATE_INDEX=true \n";


#extract chr1 region
if($chr ne ""){
#$header .= "samtools view -bh $input_bam_raw_withoutheader_super $chr > $input_bam_raw_withoutheader \n";
#$header .= "samtools index $input_bam_raw_withoutheader \n";
}


#add regroup
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/picard/AddOrReplaceReadGroups.jar I=$input_bam_raw_withoutheader O=$input_bam_raw ";
$header .= "ID=$readGroup_name LB=$readGroup_name PL=illumina PU=run SM=$sample CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT SORT_ORDER=coordinate ";
$header .= "TMP_DIR=$tmp_dir \n";

#realign target 
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/BisSNP/$BisSNPVersion -R $ref ";	
$header .= "-I $input_bam_raw ";
if($chr eq ""){
	$header .= "-T BisulfiteRealignerTargetCreator ";
}
else{
	$header .= "-T BisulfiteRealignerTargetCreator -L $chr ";
}

$header .= "-known $indel_1 ";
if($genome_version ne "mm9"){
	$header .= "-known $indel_2 ";
}

$header .= "-o $indel_target_interval -nt $cpu_cores\n ";
										
#realign 
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/BisSNP/$BisSNPVersion -R $ref ";	
$header .= "-I $input_bam_raw ";
$header .= "-T BisulfiteIndelRealigner -targetIntervals $indel_target_interval ";
$header .= "-known $indel_1 ";
if($genome_version ne "mm9"){
	$header .= "-known $indel_2 ";
}
$header .= "-compress 5 -cigar ";
$header .= "-o $input_bam_realign\n ";
										
#mark duplicate
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/picard/MarkDuplicates.jar I=$input_bam_realign O=$input_bam_realign_mdups ";
$header .= "CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=3000000 ";
$header .= "METRICS_FILE=$metrics_file ";
$header .= "TMP_DIR=$tmp_dir \n";

#baseQ recalibration
##1 countCovariant
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/BisSNP/$BisSNPVersion -R $ref ";	
$header .= "-I $input_bam_realign_mdups ";
$header .= "-T BisulfiteCountCovariates -nt $cpu_cores ";
$header .= "-knownSites $dbsnp ";
$header .= "-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate ";
$header .= "-recalFile $recalFile_before \n ";

##2 TableRecalibration
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/BisSNP/$BisSNPVersion -R $ref ";	
$header .= "-I $input_bam_realign_mdups ";
$header .= "-T BisulfiteTableRecalibration ";
$header .= "-o $input_bam_realign_mdups_recal ";
$header .= "-recalFile $recalFile_before -maxQ 40 \n ";

##3 countCovariantAfterRecalibrate
$header .= "java -Xmx$mem -jar /home/uec-00/yapingli/software/BisSNP/$BisSNPVersion -R $ref ";	
$header .= "-I $input_bam_realign_mdups_recal ";
$header .= "-T BisulfiteCountCovariates -nt $cpu_cores ";
$header .= "-knownSites $dbsnp ";
$header .= "-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate ";
$header .= "-recalFile $recalFile_after \n ";


#submit job
my $count=1;
my $outfileName = "submitPBS" . $count;
open(OUT, ">$outfileName") or die "can not open file:$!";
print OUT $header;
close OUT;
system("qsub $outfileName");
#print "$header\n";
$count++;
print "$progName is submitted\n";
system("rm $outfileName");